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ABSTRACT 


The  suitability  of  cyclostationary  spectral  analysis  for  the  detection  and 
identification  of  typical  Satcom  signals  is  examined.  An  overview  of  the  general  equations 
for  cyclostationary  spectral  analysis  is  presented.  The  cyclic  spectrum  is  derived  for  some 
common  Satcom  modulation  schemes.  Cyclic  spectrum  estimation  techniques  are 
described  with  particular  emphasis  on  the  time-smoothed  approach.  Key  aspects  of  a 
computationally-efFicient  algorithm  for  computing  the  cyclic  spectrum  estimates,  called  the 
“Fast  Fourier  Transform  (FFT)  Accumulation  Method”  (FAM)  are  outlined.  The  FAM 
algorithm  is  implemented  on  a  Digital  Signal  Processor/Personal  Computer  (DSP/PC) 
board  containing  a  TMS320c51  processor  chip  and  used  to  compute  the  cyclic  spectrum 
for  various  phase-shift-keyed  (PSK)  and  spread  spectrum  waveforms.  The  results  are 
presented  and  discussed. 


RESUME 


Ce  rapport  examine  la  pertinence  de  I’analyse  du  spectre  cyclostationnaire  pour  la 
detection  et  1’ identification  des  signaux  typiques  dans  la  communication  par  satellite.  Les 
equations  generales  pour  I’analyse  du  spectre  cyclostationnaire  sont  presentees  en  resume. 
Le  spectre  cyclique  est  derive  pour  quelques  exemples  des  modulations  utilisees  dans  la 
communication  par  satellite.  Une  description  des  techniques  de  I’estimation  du  spectre 
cyclique  est  foumie  en  soulignant  la  methode  de  la  moyenne  du  temps.  Les  aspects 
importants  d’un  algorithme  efficace  en  nombre  d’operations  pour  1’ estimation  du  spectre 
cyclique  sont  presentes.  L’algorithme  s’appelle  “La  Methode  de  I’Accumulation  des 
Transformes  de  Fourier  Rapides”  (FFT  Accumulation  Method).  L’algorithme  est  mise  en 
oeuvre  sur  un  unite  de  traitement  des  signaux  numeriques,  TMS320c51.  Le  spectre 
cyclique  est  evalue  pour  des  signaux  transmis  par  deplacement  de  phase  et  les  signaux 
etales  de  spectre.  Les  resultats  sont  presentes. 
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EXECUTIVE  SUMMARY 


With  further  increases  in  the  number  of  satellites  and  in  the  deployment  of  both 
large  and  small  terminals  expected  in  order  to  meet  Satcom  requirements  for  military 
communications,  the  risk  of  electronic  interference  causing  Satcom  stress*  is  of  concern. 
This  has  lead  to  a  desire  to  develop  effective  techniques  to  monitor  Satcom  links  for 
electronic  interference.  Cyclostationary  spectral  analysis  is  examined  to  determine  its 
suitability  for  detection  and  identification  of  typical  Satcom  signals. 

The  underlying  principle  of  cyclostationary  spectral  analysis  involves  treating 
signals  under  examination  as  cyclostationary  processes.  For  signals  with  periodic  features, 
the  advantage  of  using  a  cyclostationary  model  is  that  non-zero  correlation  occurs 
between  certain  frequency  components  of  the  signal  where  the  frequency  separation  is 
related  to  the  periodicity  of  interest.  The  general  equations  for  cyclostationary  spectral 
analysis  are  the  cyclic  autocorrelation  and  the  cyclic  spectral  density. 

Estimates  of  the  cyclic  spectrum  using  a  time- smoothed  approach  can  be  obtained 
as  follows:  frequency  components  for  a  signal  x(t)  are  evaluated  over  a  time  window,  T^, 
along  the  entire  observation  interval.  At,  at  each  window,  two  components  are  multiplied; 
and  the  resulting  sequence  of  products  is  integrated  over  time'.  This  process,  however, 
requires  significant  computational  resources.  A  computationally-efficient  algorithm  for 
cyclic  spectrum  estimation  is  called  the  “Fast  Fourier  Transform  (FFT)  Accumulation 
Method”  (FAM).  The  FAM  algorithm  computes  cyclic  spectrum  estimates  a  block  at  a 
time  by  performing  an  FFT  on  the  product  sequence  rather  than  an  integration. 

Decimation  and  data  tapering  windowing  considerations  also  reduce  the  number  of 
computations. 

The  FAM  algorithm  is  implemented  on  a  Digital  Signal  Processor/Personal 
Computer  (DSP/PC)  board  containing  a  TMS320c51  processor  chip  and  used  to  compute 
the  cyclic  spectrum  for  binary-phase-shift-keyed  (BPSK),  quadrature-phase-shift-keyed 
(QPSK),  frequency-hopped/M-ary  frequency-shift-keyed  (FH/MFSK),  and 
FH/phase-shift-keyed  (FH/PSK)  waveforms.  Periodic  components  such  as  the  carrier 
frequency,  symbol  rate  and  hop  rate  appear  as  distinct  cycle  frequency  components  in  the 
cyclic  spectrum.  Results  also  demonstrate  that  two  or  more  signals  can  be  detected  in  the 
same  band  of  interest.  Further  improvements  in  estimate  reliability  for  this  implementation 
can  be  achieved  by  considering  longer  observation  times  (i.e.  collecting  more  data 
samples).  In  turn,  strategies  including  parallel  implementations,  algorithm  development, 
and  making  use  of  advances  in  computing  power  will  help  to  minimize  compromises  in 
computation  time. 


’Stress  refers  to  any  reduction  in  communications  capability  [1] 
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1.0  INTRODUCTION 


Satellite  communication  (Satcom)  requirements  for  military  communications  have 
increased  significantly  in  recent  years  due  to  a  desire  to  increase  and  maximize  system 
throughput  and  to  implement  higher  bandwidth  and  mobile  applications.  Advances  in 
these  areas  have  led  to  an  increase  in  the  number  of  satellites  as  well  as  the  deployment  of 
both  large  and  small  Satcom  terminals.  A  potential  consequence  of  this  trend  is  the 
increased  risk  of  electronic  interference.  The  effects  of  electronic  interference  on  Satcom 
links  include  increased  noise  in  the  payload  receiver  or  ground  terminal,  a  reduction  in 
available  transponder  power  output,  and  may  ultimately  deny  timely  authorized  accesses. 

A  recent  “The  Technical  Cooperation  Program”  (TTCP)  STP-6  workshop  report  [1] 
acknowledged  electronic  interference  as  a  major  source  of  Satcom  stress’  and  looked  at 
potential  solutions.  Ideally,  solutions  which  not  only  detect  the  interference,  but  isolate 
the  interference  from  valid  Satcom  signals  are  preferred. 

Since  the  development  of  the  theoretical  concepts  by  Gardner  [2,3]  in  the  early-to- 
mid  80s,  extensive  work  has  been  carried  out  to  explore  the  potential  of  cyclostationary 
spectral  analysis  for  traditional  signal  processing  tasks.  Studies  in  potential 
cyclostationary  spectral  analysis  applications  include  time-difference-of-arrival  (TDOA) 
estimation  [4],  signal  detection  and  identification  [5-7],  and  parameter  estimation  [7,8]. 
The  underlying  principle  of  cyclostationary  spectral  analysis  involves  treating  the  signal 
under  examination  as  a  cyclostationary  rather  than  stationary  process.  For  signals  with 
periodic  features  (e.g.  digital  communications  signals),  the  advantage  of  using  a 
cyclostationary  model  is  that  non-zero  correlation  is  exhibited  between  certain  frequency 
components  when  their  frequency  separation  is  related  to  the  periodicity  of  interest. 

A  major  drawback  of  cyclostationary  spectral  analysis  techniques  is  the 
tremendous  processing  requirement.  As  a  result,  a  comparable  effort  has  been  put 
forward  to  develop  computationally-efficient  algorithms  for  computing  the  cyclic  spectrum 
estimates.  These  include  both  time-  and  frequency-smoothing  approaches  [3,9,10].  Of 
particular  interest  in  this  report  is  a  time-smoothing  algorithm  called  the  “Fast  Fourier 
Transform  (FFT)  Accumulation  Method”  (FAM)  [10].  It  takes  advantage  of  the  FFT  to 
compute  the  cyclic  spectrum  estimates  a  block  at  a  time.  Decimation  and  data  taper 
windowing  considerations  also  contribute  to  the  computational  efficiency  of  the  FAM 
algorithm. 

Several  reports  relating  to  cyclostationaiy  spectral  analysis  work  at  the  Defence 
Research  Establishment  Ottawa  (DREO)  or  by  DREO  scientists  over  the  last  5  years  have 
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also  illustrated  the  potential  of  cyclostationary  spectral  analysis  for  many  applications. 
These  include  an  examination  of  cyclic  spectrum  estimate  computations  and  their 
implications  in  a  signal  search  and  detect  environment  [1 1],  an  implementation  of  a 
computationally-efficient  algorithm  for  cyclic  spectrum  estimation  [12],  and  a  report 
featuring  a  real-time  optoelectronic  implementation  of  a  new  algorithm  which  is  able  to 
detect,  separate,  and  characterize  spread  spectrum  signals  [13], 

This  report  examines  cyclostationary  spectral  analysis  as  it  is  applied  to  detection 
and  identification  of  some  common  Satcom  signals.  An  overview  of  the  theory  for 
cyclostationary  spectral  analysis  is  presented  in  Section  2.0.  Cyclic  spectral  densities  for 
phase-shift-keyed  (PSK)  and  spread  spectrum  waveforms  are  derived.  Section  3.0 
describes  the  time-smoothed  cyclic  spectrum  measurement  technique  and  introduces  key 
concepts  of  the  FAM  algorithm.  Principal  components  of  the  experimental  setup  are 
identified  in  Section  4.0,  including  the  FAM  implementation  on  a  TMS320c5x  digital 
signal  processor  (DSP).  Results  of  the  FAM  computations  are  presented  in  Section  5.0 
for  PSK,  quadrature-phase-shift-keyed  (QPSK),  frequency-hopped/M-ary  frequency-shift 
keyed  (FH/MFSK),  and  frequency-hopped/phase-shift-keyed  (FH/PSK)  waveforms. 
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2.0  CYCLOSTATIONARY  SPECTRAL  ANALYSIS 
THEORY  OVERVIEW 


Cyclostationary  spectral  analysis  differs  from  conventional  spectral  analysis 
because,  as  the  name  implies,  signals  under  examination  are  modelled  as  cyclostationary 
rather  than  stationary  processes  [3,8],  For  the  analysis  of  waveforms  with  some  periodic 
feature(s),  a  key  advantage  of  using  a  cyclostationary  model  is  that  non-zero  correlation  is 
exhibited  between  certain  frequency  components  when  their  frequency  separation  is 
related  to  the  periodicity  of  interest.  The  value  of  that  frequency  separation  is  referred  to 
as  the  cycle  frequency.  The  non-zero  correlation  characteristic  is  reflected  in  the  cyclic 
autocorrelation  function,  denoted,  and  its  frequency-domain  equivalent,  the 

cyclic  spectral  density  or  cyclic  spectrum,  [3,8, 14].  These  are  the  general 
equations  for  cyclostationary  or  cyclic  spectral  analysis.  They  are  defined 

1 

Km  *<' ^  <i) 

-I 

2 

and 

s;(/)  4 

respectively,  where  a  =  cycle  frequency,  ^{•}=  Fourier  transform  operator,  and 

X^{f)t  J  x(u)e'^^^‘‘du 
■^'2 

The  additional  cycle  frequency  variable  in  (2)  is  reflected  in  a  two-dimensional 
representation  of  S^if)  referred  to  as  the  (f,a)  or  bifrequency  plane  (Figure  2. 1(a)).  For 
a  bandlimited  signal,  the  range  of  values  of  f  and  a  for  which  5“(/)  exists  is  referred  to  as 
the  region  of  support.  Figure  2. 1(b)  illustrates  the  region  of  support  for  lowpass  signals. 
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(a)  (b) 


Figure  2.1  (a)  Bifrequency  plane  (b) 
Region  of  support  for  lowpass  signal 


For  waveforms  which  have  undergone  common  signal  processing  operations  (e.g. 
filtering,  modulation),  evaluation  of  (1)  and  (2)  can  be  cumbersome.  However,  in  [3,8] 
several  properties  or  relationships  are  described  which  help  simplify  this  task.  One  useful 
result  is  the  expression  of  the  cyclic  autocorrelation  for  the  product  of  two  independent 
waveforms.  It  is  shown  [3,8],  that  for  x(t)=d(t)c(t), 

«“(V=E  (3) 

P 


Thus,  the  cyclic  autocorrelation  of  a  product  of  two  independent  waveforms  is  given  by 
the  discrete  convolution  in  cycle  frequency  of  the  independent  cyclic  autocorrelation 
functions.  The  cyclic  spectrum  is  then  given  by 

-0&  ^ 

Note  that  the  cyclic  spectrum  is  continuous  in  f  and  discrete  in  a.  The  results  of  (3)  and 
(4)  are  used  to  evaluate  the  cyclic  spectrum  of  some  typical  Satcom  signals. 

2.1  Cyclostationary  Spectral  Analysis  of  Satcom  Signals 

Common  Satcom  signals  may  include  PSK  and  QPSK  modulations  depending  on 
data  rate  requirements.  FSK  modulation  schemes  may  be  used  for  low  data  rate 
applications  and  in  cases  where  bandwidth  efficiency  is  of  less  importance.  In  addition, 
direct  sequence  code  division  multiple  access  (CDMA),  FH/MFSK  and  frequency-hopped 
DPSK  (FH/DPSK)  are  now  becoming  more  commonplace  for  a  number  of  reasons,  not 
the  least  of  which  is  their  robustness  to  electronic  interference.  For  detection  and 
identification  applications,  it  is  thus  appropriate  to  individually  examine  the  cyclic 
spectrum  of  these  more  typical  Satcom  signals. 
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2.1.1  PSK 


Let  the  PSK  signal  be  described  by 

x{t)  =^d{i)c{t) =d{t)  cosilTifJ  +  <1)^) 

where  d(t)  is  the  data  signal  and  c(t)  is  the  cosine  or  carrier  term.  Since  d(t)  and  c(t)  are 
independent  time  series,  the  cyclic  autocorrelation  o^x(^)  may  be  given  by  (3).  Applying 
the  cyclic  autocorrelation  function  definition,  (1),  to  the  cosine  term  c(t)  yields 

a=0 

«=±2/, 

othemise 

The  cyclic  spectrum  for  the  cosine  term  is  hence  given  by 

i'/C/)'  a-±2/<, 

0 ;  otherwise 

which  consists  of  delta  functions  at  the  carrier  frequency  for  a=0  and  at  /=0  for 
a=±lf^.  The  cyclic  spectrum  fox  x(t)  is  then  simply  the  cyclic  spectrum  of  the  data 
sequence,  d(t),  replicated  at  the  cosine  term  delta  functions. 

a=±2/„ 


‘/2COS(27t4t); 

.0; 


Now  consider  the  data  signal  d(t)  which  may  also  be  rewritten  as  a  time-sampled 
random  data  waveform,  a(t),  with  the  time  sampling  described  by  a  periodic  data  carrier 
XQVxxi,p(t),  of  period 

dit)=a(t)p(t)=  ainTj6it-nT^)=±  ^ 

n=-oc 


Because  a(t)  and  p(t)  are  statistically  independent  we  may  again  obtain  the  cyclic 
autocorrelation  function  from  (3).  Assuming  a(t)  is  purely  stationary,  noting  that  the 
Fourier  transform  of  p(t)  consists  of  delta  functions  at  harmonics  of  the  symbol  rate  with 
amplitude  and  performing  considerable  algebraic  manipulation  yields, 


^otherwise 
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where  k  and  /  are  integers.  It  is  observed  that  d(t)  is  a  basic  component  for  many  digital 
modulation  formats,  which  suggests  that  the  cyclic  spectrum  of  d(t)  will  appear  in  a 
number  of  other  digital  modulation  examples. 

2.1.2  QPSK 

A  QPSK  waveform  may  be  described  by 

x{t)=d,{t)  cos{2ufj) -dg  (0  sin  (27r//) 

with  the  data  signals  d/t)  and  dgO)  as  described  above  for  d(t).  The  derivation  of  the 
cyclic  spectrum  for  the  QPSK  signals  follows  a  similar  approach  as  for  the  PSK  signal. 
Unlike  the  PSK  example  however,  the  cyclic  autocorrelation  for  QPSK  yields  non-zero 
correlation  only  at  the  carrier  frequency.  The  cyclic  spectrum  for  a=±2fg  tends  to  zero 
when  dj(t)  and  dgft)  are  balanced  -  the  quadrature  components  cancel.  Letting 
dft)=dg(t)=d{t),  the  QPSK  cyclic  spectrum  may  be  described  by 

where  Sf(f)  is  as  given  for  PSK. 

2.1.3  FH/MFSK 


From  [14],  the  cyclic  spectrum  for  multiple  FH/MFSK  hops  can  be  obtained  by 
first  considering  the  case  of  a  single  hop, 

x(t)=d{l)e 

where  d(t)  is  as  described  for  the  PSK  example  except  with  a(nTJ=0,J .  The  cyclic 
spectrum  is 


0; 


a=0 


othen\’ise 


where  Ti  =  hop  period,  Sinc(x)  =  sin( 7dc)/( nx).  The  envelope  of  the  cyclic  spectrum  is  a 
Sinc(x)  fiinction  in  a  at  the  carrier  frequency.  For  the  case  of  multiple  hops,  a  similar 
envelope  is  expected  at  each  of  the  hop  carrier  frequencies  along  with  a  component  at  a 
cycle  frequency  equal  to  the  difference  in  the  hop  carrier  frequencies, 
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s;{f)- 


—  T^ 
At  " 

1  -.2 


At 


T^Sinc^ 


«=±(/;-/j) 


6/=/j-/2 


0; 


otherwise 


However,  the  cyclic  autocorrelation  for  any  given  hop  or  between  any  number  of  hops  will 
be  non-zero  for  a  given  (/"^  arj  only  for  the  duration  of  a  hop.  The  cyclic  spectrum  thus 
tends  to  zero  for  increasing  time  window,  At. 

2.1.4  FH/PSK 

The  FH/PSK  waveform  is  described  by 

x{t)-m{i)d(t)e 

where  m(t)  =  ±1  is  the  PSK  data  and  d(t)=  0,1  is  the  hop  on/oflf  pulse  as  given  above. 

From  [14],  the  spectrum  is  described  by 

where  A'=  Fourier  coefficients  of  x(t)  evaluated  over  a  hop  period,  (i.e.  a  time-limited 
BPSK  signal).  S“(f )  is  only  non-zero  for  cycle  frequencies  equal  to  integer  multiples  of 
the  hop  rate.  As  with  the  FH/MFSK  case,  the  cyclic  spectrum  tends  to  zero  for  increasing 
At. 
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3.0  CYCLIC  SPECTRAL  DENSITY  ESTIMATION 


In  practice,  signals  under  consideration  are  defined  over  a  finite  time  interval.  As  a 
result,  the  cyclic  spectral  density  cannot  be  measured  exactly.  Estimates  of  the  cyclic 
spectral  density  can  be  obtained  via  time-  or  frequency-smoothing  techniques  [3,8-10].  In 
this  report,  a  time-smoothing  approach  is  investigated  where  the  estimate  of  the  cyclic 
spectral  density  is  given  by  the  time-smoothed  cyclic  periodogram  [3], 
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where  At=obser\>ation  time,  T^=time  window  length,  and 

Xj.  (t,/)^  short  time  Fourier  transform  of  x{t) 
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Figure  3.1  illustrates  this  process  graphically  for  an  arbitrary  (f^  a„).  For  some  signal,  x(t), 
frequency  components  are  evaluated  over  a  time  window,  along  the  entire  observation 
interval  At.  At  each  window,  two  components,  centred  about  some/,  and  separated  by 
some  are  multiplied  and  the  resulting  sequence  of  products  is  integrated  over  time. 


Spectral  components  obtained  from  the  short-time  Fourier  transform  have  a 
resolution,  Af  of  approximately  V/  ,.  The  product  of  two  such  components  also 
provides  a  resolution  capability  in /of  for  the  cyclic  spectrum  estimates.  To 
minimize  the  effects  of  spectral  and  cyclic  leakage  introduced  through  frequency 
component  sidelobes  [10,14],  a  data  taper  window  (e.g.  Hamming  window)  is  applied. 
Estimate  reliability  is  further  enhanced  if  a  substantial  amount  of  smoothing  is  carried  out 
over  At  to  reduce  the  variance  of  the  estimate.  From  Grenander’s  Uncertainty  Condition 
[3],  this  requires  that 
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Figure  3.1  Cyclic  spectrum  estimation  via  time-smoothing  (a) 
Short-time  Fourier  transform  over  At  (b)  Sequence  of 
frequency  products  from  short-time  Fourier  transforms 


which  means  that  the  observation  time  must  greatly  exceed  the  time  window  used  to 
compute  the  spectral  components.  Furthermore,  it  is  shown  [3]  that  the  cycle  frequency 
resolution  of  the  estimate  is  Aa  »  This  implies  that  an  estimate  for  some  (f„  cc^ 

represents  a  very  small  area  on  the  bifrequency  plane  (Figure  3.2)  and  to  adequately 
represent  the  cyclic  spectrum  requires  a  significant  number  of  estimates.  It  follows  that 
using  a  brute  force  approach  to  obtain  estimates  is  computationally  demanding. 


Figure  3.2  One  cyclic  spectrum  estimate  on  a 
biffequency  plane  for  a  lowpass  signal 


3.1  The  FFT  Accumulation  Method 

Computationally-efFicient  techniques  for  obtaining  cyclic  spectrum  estimates  have 
been  developed  to  try  to  reduce  the  number  of  computations  and  are  described  in 
[9,10,14].  One  particular  algorithm  based  on  the  time-smoothing  approach  is  called  the 
FFT  Accumulation  Method  (FAM).  The  strategy  of  the  FAM  algorithm  is  to  divide  the 
bifrequency  plane  into  smaller  areas  called  channel  pair  regions  and  compute  the  estimates 
a  block  at  a  time  using  the  efficient  FFT. 

The  first  step  in  describing  the  FAM  algorithm  is  to  express  (5)  in  discrete  terms 
yielding. 


where 
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and  M>[n]=  data  taper  window.  The  discrete  equivalents  of/ and  a  are  k  and  y 
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respectively.  Parameters  N  and  N’  represent  the  length  of  the  sequence  and  length  of  the 
discrete  short-time  FFT.  Extending  Grenander’s  Uncertainty  Condition,  the  reliability  of 
the  discrete  estimate  now  depends  on, 

In  the  FAM  algorithm,  spectral  components  of  a  sequence,  xfnj,  are  computed 
using  (8).  Two  components  are  multiplied  to  provide  a  sample  of  a  cyclic  spectrum 
estimate  representing  a  finite  area  on  the  bifrequency  plane  called  a  channel  pair  region 
(Figure  3.3(a)).  A  sequence  of  samples  for  that  particular  area  can  be  obtained  by 
multiplying  the  same  two  components  for  a  series  of  consecutive  short-time  FFTs  which 
span  the  length  of  the  input  sequence.  However,  instead  of  computing  an  average  of  the 
product  sequence,  as  in  (7),  an  FFT  is  performed  on  the  second  sequence  to  further 
resolve  the  estimates  along  the  cycle  frequency  axis.  The  second  FFT  generates  a  block  of 
cyclic  spectrum  estimates;  one  of  these  (i.e.  the  DC  component)  would  include  the 
estimate  obtained  by  time-averaging  (Figure  3.3(b)).  Estimates  for  the  entire  cyclic 
spectrum  are  obtained  by  following  this  procedure  for  all  combinations  of  spectral 
component  pairs  from  the  short-time  FFTs  (i.e.  channel  pair  regions)  as  shown  in  Figure 
3.3(c). 


(a)  (b)  (c) 

Figure  3.3  Partitioning  of  the  bifrequency  plane  (a)  Channel  pair  region  represented 
by  the  product  of  two  spectral  components  (b)  Block  of  cyclic  spectrum  estimates 
for  one  channel  pair  region  from  the  FFT  of  the  product  sequence  (c)  Channel  pair 
regions  of  the  entire  bifrequency  plane 


From  Figure  3.3(b),  it  is  obvious  that  estimates  toward  the  “top”  and  “bottom”  of 
the  channel  pair  region  do  not  satisfy  Grenander’s  Uncertainty  Condition.  To  address  this 
problem,  one  approach  involves  taking  only  those  resolved  cyclic  spectrum  components 
that  are  the  equivalent  in  y  of  ±^r«/2  away  from  the  centre  of  the  channel  pair  region 
[10].  However,  taking  only  some  of  the  cyclic  spectrum  estimates  will  leave  gaps  in  the 
bifrequency  plane  (Figure  3.4).  This  can  be  resolved  by  exploiting  the  effect  of  applying  a 
data  taper  window  in  (8)  on  the  frequency  resolution  to  obtain  better  channel  pair  region 
coverage  as  shown  in  Figure  3.5. 


Figure  3.4  Gaps  in  bifrequency  plane  when 
considering  only  estimates  that  are  ±KeJ2 
away  from  the  channel  pair  region  centre 


By  computing  a  block  of  estimates  using  an  efficient  FFT,  the  number  of 
computations  is  significantly  reduced.  Additional  savings  in  computational  requirements 
can  be  achieved  by  reducing  the  number  of  short-time  FFTs.  This  is  accomplished  by 
decimating  the  product  sequence  and  performing  short-time  FFTs  every  Z  samples  instead 
of  at  every  sample.  The  parameter  L  is  referred  to  as  the  overlap  parameter.  It  is  shown 
[10,1 1,14]  that  to  guard  against  aliasing  and  cycle  leakage  effects  on  the  cyclic  spectrum 
estimates,  /4.  A  comparison  of  the  computational  requirements  of  the  FAM 
algorithm  show  a  reduction  by  a  factor  of  over  a  brute  force  approach  [10,14]. 
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Taking  into  consideration  the  issues  discussed  above,  the  FAM  algorithm  can  be 
expressed  by  the  following, 

Y  +cjAy  ^  Y  Y 

s  ",  {’’Ah  =  E  X^(„+mL,k+l^)X;{n+mLA~)  ” 

where, 

Y'-l 

Xf,,/(r,k)- M>\d\x{a+r\e  ^  ,w\d\— window  coefficients 

fl=o  (10) 

and  P  =  length  of  the  product  sequence. 
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4.0  IMPLEMENTATION  AND  TEST  SETUP 


The  purpose  of  this  section  is  to  present  the  components  of  the  experimental  setup 
used  to  test  an  implementation  of  the  FAM  algorithm.  The  FAM  algorithm  was 
implemented  on  a  PC/c5x  DSP  board  made  by  Loughborough  Sound  Images  Inc.  The 
DSP  board  includes  a  Texas  Instruments  TMS320c51  DSP  chip  which  performs  the  FAM 
computations.  A  general  description  of  the  implementation  is  included. 


4.1  Experimental  Setup 


The  key  components  in  the  experimental  setup  are  as  follows:  signal  generation, 
received  signal  processing,  baseband  signal  processing  and  data  collection,  cyclic  spectrum 
estimate  computation,  and  hardware/user  interface.  They  are  illustrated  in  Figure  4.1. 


In  order  to  simplify  the  setup  without  affecting  the  evaluation  of  the  FAM 
algorithm,  input  test  signals  are  generated  at  IF  and  combined.  In  effect,  these  signals 
represent  a  typical  front  end  of  a  communications  receiver.  The  resulting  signal  is 
translated  via  quadrature  detection,  producing  the  I  and  Q  baseband  signals  to  be 
analyzed.  An  analog-to-digital  (A/D)  board  collects,  digitizes  and  stores  the  input 
waveform.  At  the  command  of  the  DSP,  the  stored  waveform  is  read  into  DSP  memory 
for  the  FAM  computations.  Further  details  on  each  one  of  these  blocks  are  described  in 
[14]. 
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4.2  FAM  implementation  description 
4.2.1  FAM  computations  overview 

The  FAM  algorithm  was  implemented  in  assembler  for  the  TMS320c5 1  DSP 
board  made  by  Loughborough  Sound  Images,  Inc.  The  assembler  program  is  divided  into 
three  modules:  data  pre-processing,  computation  of  frequency  components,  and 
computation  of  cyclic  spectrum  estimates.  A  flow  chart  showing  the  main  functions  of  the 
assembler  program  according  to  these  three  modules  is  given  in  Figure  4.2.  In  this  report, 
a  summary  of  the  computational  process  is  presented.  Specific  details,  however,  are 
contained  in  [14]. 

At  the  start  of  the  assembler  program  for  the  FAM  computations,  parameters 
including  the  input  sequence  length  (A'^=1024),  short-time  Fourier  transform  length 
(N-32),  and  overlap  parameter  (Z  =  8)  are  initialized.  A  master  reset  of  the  A/D  board  is 
initiated  by  the  DSP.  When  the  DSP  receives  a  ready  signal  from  the  A/D  board,  the  data 
collected  on  the  A/D  board  is  read  into  DSP  board  memory.  A/D  samples  are  represented 
as  12-bit  signed  numbers  using  offset  binary  coding  and  need  to  be  converted  to  2's 
complement  representation  for  the  DSP  computations.  In  this  implementation,  the  input 
sequence  is  zeropadded  so  that  the  length  of  the  second  sequence,  P,  will  be  a  power-of-2 
to  facilitate  the  use  of  a  Radix-2  FFT  routine. 

A  series  of  overlapping  short-time  FFTs  is  performed  to  obtain  the  spectral 
components  of  the  input  sequence.  The  output  of  each  short-time  FFT  is  stored  in  an 
array.  All  frequency  components  are  normalized  to  16  bits  to  reduce  the  memory 
requirements  when  products  are  subsequently  formed  from  these  frequency  components 
and  stored. 

Complex  product  sequences  for  a  particular  channel  pair  region  at  (k^  yj  are 
formed  by  multiplying  the  appropriate  spectral  components  together  over  all  the  short- 
time  FFTs  or  time  hops.  For  each  product,  a  time-shift  exponential  factor  is  applied  to 
preserve  the  coherence  of  the  spectral  components  between  successive  short-time  FFT 
hops.  A  second  FFT  is  performed  on  the  product  sequence  to  obtain  the  cyclic  spectrum 
estimates  for  that  channel  pair  region.  As  described  in  Section  3.1,  only  estimates 
representing  the  area  res  12  away  from  the  centre  of  the  channel  pair  region  are  retained 

to  maintain  statistical  reliability.  These  estimates  are  then  stored  in  DSP  memory.  This 
process  is  repeated  for  all  possible  combinations  of  spectral  components  or  channel  pair 
regions.  The  DSP  then  sends  a  signal  to  the  PC  to  indicate  that  the  cyclic  spectrum 
estimates  are  ready  to  be  transferred.  It  is  noted  [10]  that  since  the  product  sequences  and 
the  second  FFT  calculations  in  the  third  module  are  independent  of  one  another,  they  can 
be  implemented  in  parallel,  thus  saving  some  computation  time. 
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4.2.2  DSP/PC  Handshaking 

The  DSP  board  is  contained  in  an  IBM-compatible  PC  and  interfaces  to  the  PC 
through  I/O  and  memory-mapped  (status  and  control)  registers.  The  PC  loads  the  DSP 
program  and  resets  the  DSP  board  to  start  the  FAM  computations.  When  the  DSP 
indicates  that  the  computations  have  been  completed,  the  PC  transfers  and  formats  the 
results  of  the  FAM  computations,  and  produces  an  ASCII  file  of  the  formatted  output. 
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5.0  EXPERIMENTAL  RESULTS 


Using  the  FAM  implementation  described  in  the  previous  sections,  the  cyclic 
spectrum  was  estimated  for  different  cases  of  experimental  and  simulated  input  data.  The 
results  are  presented  for  some  common  modulation  schemes  including  BPSK,  QPSK, 
FH/MFSK,  and  FH/PSK. 

To  verify  the  implementation  of  the  FAM  algorithm,  two  CW  signals,  with  carrier 
frequencies  close  to  one  another,  were  generated  as  the  input  signal.  Figure  5.1  illustrates 
the  resulting  cyclic  spectrum.  A  cyclic  component  is  observed  at  cc^  which  corresponds  to 
the  difference  between  the  two  carrier  frequencies.  This  illustrates  the  potential  for  cyclic 
spectral  analysis  to  be  used  to  detect  two  or  more  signals  within  the  same  band  of  interest 
where  it  would  be  difficult  using  the  conventional  power  density  spectrum  (along  a  =0). 


a 


(a)  (b) 

Figure  5.1  Cyclic  spectrum  for  two  CW  signals  (a)  CW  signals  in  the  same  frequency 
bin  (b)  Same  frequency  bin  -  side  view 


For  BPSK  waveforms,  a  comparison  of  the  cyclic  spectra  for  low  and  high  data 
rate  signals  in  AWGN  is  given  in  Figure  5.2.  From  the  PSK  example  of  Section  2. 1 . 1,  a 
cyclic  component  at  the  data  is  expected.  The  low  data  rate  BPSK  signal,  however, 
shows  no  distinguishable  component  in  the  cyclic  spectrum.  With  the  relatively  short 
input  sequence  length  implemented  (A^=1024)  compared  to  the  symbol  period,  there  was 
not  enough  signal  captured  to  obtain  an  accurate  estimate  for  small  values  of  a.  This 
result  illustrates  that  using  as  large  an  observation  time  window  as  is  practical  would  be 
preferred  for  detection  applications.  In  the  case  of  the  higher  data  rate  signal,  a  cyclic 
component  is  visible  as  more  data  symbols  are  included  in  the  calculation  for  the  same 
observation  time  owing  to  the  shorter  bit  period. 
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(C)  (f) 

Figure  5.2  Cyclic  spectmm  for  a  BPSK  signal  in  AWGN  (a)  E^/Ng  >  30dB,  /fj  = 
128kbps  (b)  E^/Ng  =  lOdB,  =  128kbps  (c)  E^/Ng  =  OdB,  =  128kbps  (d)  E^/Ng  > 
30dB,  Rt  =  1 .024  Mbps  (e)  E^  /Ng  =  lOdB,  =  1 .024  Mbps  (f)  £*  /Ng  =  OdB,  Rt  = 
1.024  Mbps 


The  cyclic  spectrum  for  a  QPSK  waveform  was  computed  and  is  shown  in  Figure 
5.3.  The  symbol  rate  of  the  signal  is  exhibited  in  the  bifrequency  plane.  However,  no 
spectral  correlation  associated  with  the  carrier  frequency  is  observed  as  in  the  case  of  the 
BPSK  waveform  (Figure  5.3(a)).  As  noted  in  Section  2.1.2,  this  component  is  cancelled 
out  by  the  in-phase  and  quadrature  components  of  the  QPSK  waveform  which  allows  the 
two  modulation  schemes  to  be  distinguished  from  one  another.  The  non-zero  values 
surrounding  the  cyclic  component  are  related  to  the  variance  of  the  estimate  for  the 
relatively  short  input  sequence. 


Figure  5.3  Cyclic  spectrum  for  BPSK  and  QPSK  waveforms  (a)  Real  BPSK  signal, 
Rf,= 1.024Mbps  (b)  Real  QPSK  signal,  R=512ksps 


Figure  5.4  shows  the  cyclic  spectrum  computed  for  an  FH/MFSK  signal.  There 
are  several  hops  observed  in  the  frequency  band  of  interest  in  the  form  of  distinct  peaks 
along  the  a  =0  axis.  However,  no  spectral  correlation  is  observed  between  the  carrier 
frequencies  of  the  hops,  as  was  suggested  in  Section  2.1.3,  because  in  this  example,  only 
one  transmission  over  a  particular  hop  period  is  received.  As  a  result,  when  the 
overlapping  short-time  FFTs  are  computed  in  the  FAM  algorithm,  two  or  more  frequency 
hops  are  unlikely  to  be  captured  simultaneously.  Consequently,  the  only  non-negligible 
products  in  the  second  sequence  are  those  centred  about  the  individual  hops.  However, 
for  multiple  transmissions  over  a  hop  period,  correlation  would  also  be  observed  at  a  cycle 
frequency  equal  to  the  difference  in  carrier  frequency. 

Similar  cyclic  spectra  are  observed  in  Figure  5.5  for  FH/PSK  and  burst  PSK 
waveforms.  Although  Figures  5.5(a)  and  5.5(c)  do  not  reveal  any  distinct  cyclic  features, 
the  profiles  show  the  periodic  component  of  the  signal,  the  hop  rate.  The  cyclic 
components  of  the  FH/PSK  example  appear  “noisier”  than  in  the  burst  PSK  example  due 
to  the  fact  that  they  are  computed  from  a  single  hop.  The  profile  also  shows  the  effect  of 
the  superposition  of  the  5/«cr-like  amplitude  response  function  associated  with  the  Fourier 
transform  of  the  hop  data. 
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(C)  (d) 


Figure  5.5  Cyclic  spectrum  for  burst  PSK  and  FH/PSK  waveforms  (a)  Burst  PSK 
waveform  (b)  Burst  PSK  waveform  -  side  view  (c)  FH/PSK  waveform  (d)  FH/PSK 
waveform  -  side  view 
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6.0  SUMMARY  AND  CONCLUSIONS 


Cyclostationary  spectral  analysis  is  a  technique  which  exploits  the  periodicities 
inherent  in  digitally-modulated  signals  and  enables  signals  with  overlapping  temporal 
spectra  to  be  isolated  at  their  cycle  frequency.  It  involves  treating  the  signal  under 
examination  as  a  cyclostationary  rather  than  stationary  process.  The  basic  equations  for 
cyclostationary  spectral  analysis  (i.e.  the  cyclic  autocorrelation  and  cyclic  spectral  density) 
are  reviewed.  Cyclic  spectra  are  evaluated  for  typical  Satcom  signals  including  PSK  and 
spread  spectrum  waveforms. 

Techniques  for  estimating  the  cyclic  spectrum  are  derived  from  time-  and 
frequency-smoothing  approaches.  Obtaining  high  reliability  estimates,  however,  requires 
significant  computational  resources.  This  report  considers  a  time-smoothing  approach  and 
examines  the  key  components  of  the  computationally-efficient  FAM  algorithm.  The 
reduction  in  the  number  of  computations  is  attributed  to  the  use  of  an  FFT  to  compute 
cyclic  spectrum  estimates  a  block  at  a  time  as  well  as  decimation  and  data  taper 
windowing  considerations.  The  FAM  algorithm  is  implemented  on  a  TMS320c5x  DSP  in 
three  main  processes:  data  pre-processing,  computation  of  frequency  components,  and 
computation  of  cyclic  spectrum  estimates.  The  DSP  is  then  used  to  compute  the  cyclic 
spectra  of  BPSK,  QPSK,  FH/MFSK,  and  FH/PSK  waveforms. 

Theoretical  and  FAM  estimates  of  the  cyclic  spectrum  have  demonstrated  the 
technique’s  ability  to  identify  typical  Satcom  waveforms.  Periodic  components  such  as  the 
carrier  frequency,  symbol  rate  and  hop  rate  are  exhibited  in  the  form  of  distinct  cycle 
frequency  components  in  the  bifrequency  plane.  Results  also  indicate  that  two  or  more 
signals  could  be  detected  within  the  same  band  of  interest.  Further  improvements  in 
estimate  reliability  for  this  implementation  can  be  achieved  by  increasing  the  observation 
time,  though  at  the  expense  of  computation  time.  However,  this  can  be  overcome  by 
considering  parallel  implementations  of  the  third  process  in  the  FAM  algorithm.  This, 
along  with  advances  in  computing  power  and  algorithm  development  would  enhance  the 
potential  of  cyclostationary  spectral  analysis  to  be  a  practical  tool  for  detection  and 
identification. 
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